%% This function produces the input needed for Figure 2 in Nucera, Sarno and Zinna "Currency Risk Premiums Redux", forthcoming on the RF
function [tabtmp, txteigc, varetmp] = f_get_eigenvalues(vhatold, betahatold, MatPrt, w, pmax);


MatPrt  = MatPrt * 12;
vahtold = vhatold;

T = size(MatPrt, 1);

% Case 0: var(f) 
avgftmp       = mean(vhatold)';
[pvaltmp]     = f_avgretNW(vhatold); 
pvalftmp      = pvaltmp';
stdftmp       = nanstd(vhatold)';
srftmp        = avgftmp./stdftmp;  
[tmp, rnktmp] = sort(avgftmp, 'descend');

% Case 1: eig beta*Vf*beta'  
covftmp     = cov(vhatold);      
betaftmp     = betahatold;
eig1tmp     =  sort(real(eig(betaftmp*covftmp*betaftmp')), 'descend');
eig1tmp     = eig1tmp(1:pmax);

% Case 2: eig beta*(Vf + (1+w)*(muf*muf'))*beta' 
muftmp      = nanmean(vhatold)';
eig2tmp     = sort(real(eig(betaftmp*(covftmp + (1+w).*(muftmp*muftmp'))*betaftmp')), 'descend');
eig2tmp     = eig2tmp(1:pmax);

% compute sige
for ntmp = 1:size(MatPrt, 2)
yytmp  = MatPrt(2:end, ntmp);
xxtmp  = [vhatold];
iFin   = find(isfinite(yytmp)==1);
yytmp  = yytmp(iFin, 1); 
xxtmp  = xxtmp(iFin, :); 
iota   = ones(size(yytmp, 1), 1);
nlag    = round((4*(size(yytmp, 1)/100))^(2/9));
output  = nwest(yytmp, xxtmp, nlag);
restmp  = output.resid';
varentmp(ntmp)  = var(restmp);
clear yytmp xxtmp iota
end

varetmp = sum(varentmp);

% replace mean in the test asset returns
MatPrtAdj = MatPrt;
for ntmp = 1:size(MatPrt, 2)
avgptmp(:, ntmp) = nanmean(MatPrtAdj(:, ntmp));     
inan      = find(isnan(MatPrtAdj(:, ntmp))==1);
MatPrtAdj(inan, ntmp) = avgptmp(:, ntmp);
end


objtmp = (1/T * (MatPrtAdj'*MatPrtAdj) + w * (avgptmp'*avgptmp));
eig3tmp   = sort(eig(objtmp), 'descend');
eig3tmp   = eig3tmp(1:pmax);


eig1stdtmp = eig1tmp./varetmp;
eig2stdtmp = eig2tmp./varetmp;
eig3stdtmp = eig3tmp/varetmp;

cheig1stdtmp = [eig1stdtmp(1:end-1)-eig1stdtmp(2:end);... 
                nan(1, size(eig1stdtmp, 2))];
            
cheig2stdtmp = [eig2stdtmp(1:end-1)-eig2stdtmp(2:end);... 
                nan(1, size(eig2stdtmp, 2))];            

cheig3stdtmp = [eig3stdtmp(1:end-1)-eig3stdtmp(2:end);... 
                nan(1, size(eig3stdtmp, 2))];            
            
            
infoT.fmt    = '%10.2f';
tabtmp       = [pvalftmp 100*avgftmp stdftmp.^2 sqrt(12)*srftmp rnktmp...
                eig1tmp       eig2tmp       eig3tmp...
                eig1stdtmp    eig2stdtmp    eig3stdtmp...
                cheig1stdtmp  cheig2stdtmp  cheig3stdtmp];

txteigr       = strcat('F_',num2str([1:pmax]'));


fprintf(['\n'])
fprintf('Figure 2: entries of Table A8 and input to Figure 2 for a given w\n')
fprintf(['\n'])
fprintf('p_eig: plain eigenvalues \n')
fprintf('n_eig: normalized eigenvalues \n')
fprintf('dn_eig: difference in normalized eigenvalues \n')
fprintf('1, 2 and 3 indicate the type(e.g., PCA, RP,..) of the sigma matrix considered (see the legend of Table A8) \n')
fprintf(['\n'])

txteigc       = strvcat('pval',    'avg(f)',  'var(f)', 'sr(f)',  'rnk(f)' ,...
                        'p_eig1',    'p_eig2',    'p_eig3', ...
                        'n_eig1',    'n_eig2',    'n_eig3',...
                        'dn_eig1',   'dn_eig2',   'dn_eig3');                    
                    
% print results
infoT.rnames = strvcat(' ', txteigr);
infoT.cnames = strvcat(txteigc);
mprint(tabtmp,infoT)